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ABSTRACT 

We test the effect of assumptions about stellar motion on the behavior of gravitational 
instabilities in protoplanetary disks around solar-type stars by performing two simulations 
that are identical in all respects except the treatment of the star In one simulation, the star 
is assumed to remain fixed at the center of the inertial reference frame. In the other, stellar 
motion is handled properly by including an indirect potential in the hydrodynamic equations 
to model the star's reference frame as one which is accelerated by star/disk interactions. The 
disks in both simulations orbit a solar mass star, initially extend from 2.3 to 40 AU with a 
vj^^''^ surface density profile, and have a total mass of 0.14 Mq. The 7 = 5/3 ideal gas is 
assumed to cool everywhere with a constant cooling time of two outer rotation periods. 

The overall behavior of the disk evolution is similar, except for weakening in various 
measures of GI activity by about at most tens of percent for the indirect potential case. Over- 
all conclusions about disk evolution in earlier papers by our group, where the star was always 
assumed to be fixed in an inertial frame, remain valid. There is no evidence for independent 
one-armed instabilities, like SLING, in either simulation. On the other hand, the stellar mo- 
tion about the system center of mass (COM) in the simulation with the indirect potential is 
substantial, up to 0.25 AU during the burst phase, as GIs initiate, and averaging about 0.9 
AU during the asymptotic phase, when the GIs reach an overall balance of heating and cool- 
ing. These motions appear to be a stellar response to nonlinear interactions between discrete 
global spiral modes in both the burst and asymptotic phases of the evolution, and the star's 
orbital motion about the COM reflects the orbit periods of disk material near the corotation 
radii of the dominant spiral waves. This motion is, in principle, large enough to be observable 
and could be confused with stellar wobble due to the presence of one or more super- Jupiter 
mass protoplanets orbiting at lO's AU. We disc uss why the excursions in our simulation are 
so much larger than those seen in simulations bv lRice et alj (l2003ah . 
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1 INTRODUCTION 



As reviewed in IPurisen et al.l <2007h . gravitational instabilities 
(GIs) will lead to the dynamic growth of spiral waves in proto- 
planetary disks whenever the Toomre parameter Q = CsK/irCE 
becomes less than about 1.5 to 1.7 or so. Here, Cs is the sound 
speed, K is the epicyclic frequency (~ SI, the disk orbital angu- 
lar frequency, in a nearly Keplerian disk), and E is the surface 
mass density. As the spirals grow to nonlinear amplitude, either 
they will saturate at a point where radiative cooling balances the 
heat generated by the waves, or they will fragment into bound 
clumps if the radiative cooling time tcooi is compar able to or shorter 
than the disk rotation period Prot ~ 2-k/Q, (e.g.. lGammiell200ll : 
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iRice et al.|[2003bl , |2005| ; JMeiia et alj|200^ : Iciarke et ai]|2007h . In 
nonfragmenting disks, mass transport rates for realistic conditions 
can result in effective Shakura & Sunyaev ( 1973) of- v iscosity pa- 
ramet ers of order 10"^ (e.g., iLodato & Ricd |2004 iBolev et al] 
}20() g), sufficient to dom inate major phases of disk evolution (e.g., 
Voro bvov & Basij|2007h . Fragmentation of disks due to fast cool- 
ing may be an important planet formation mechanism under some 
cond itions (e.g.. Boss 1997, 2000. 2001; Mayer et al. 2002. 2004|; 
Sta matellos&Whitworthl2009l ; lBolevll2009l : IClarkel2009l : lRafikovl 
120091) . 

Early simulations of GIs in disks by the Indiana Univeristy 
hydr odynamics group ( l UH G) included a hydrodynamically active 
star jPickett et al]|l998l. [20bo). but the disks we could model with 
active stars were only ten times larger than the stellar radius. To pro- 
duce more realistically sized protoplanetary disks, we introduced a 
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central hole and fixed the star at the center of our cylindrical com- 
putational grid JPickett et alj|2003l; iMeiia et al.ll2005l: iBolev et alj 



iLJB 

enl2 



l2006Ll2007l : ICail2006l . lCai et alj2008l : rBolev & Durisel^200 8^. The 
hydrodynamics equations remained written as if the grid reference 
frame were inertial. With an artificially fixed central stellar poten- 
tial in the inertial frame, star/disk interactions that might shift the 
star off the system center of mass (COM) cannot be properly mod- 
eled. Acceleration of the star and subsequent feedback into the spi- 
ral structure is explicitly suppressed. In a sense, the star has infinite 
inertia and becomes a sink for momentum. Under these co nditions, 
SLING amplification JAdams et al.ll989l ; IShu et al.ll990l) of spiral 
structure cannot be properly computed. The concern about our ear- 
lier calculations is then two-fold: 1) The physics in our simulations 
of Gl-active disks may be inaccurate, e.g., leading to erroneous es- 
timates of the effective as for mass transport, and may lack impor- 
tant features, e.g., stronger and more coherent one-armed spirals. 
2.) The stellar motions that we suppress could be a significant and 
observable signature of GI activity in disks (Rice et al. 2003a). 

Authors have approached the problem of stellar motion 
for disk simulations in different ways. In Smoothed Particle 
Hydrodynamics (SPH) simulations of GIs (Rice et al. 2003b; 
Mayer et al. 2004; Lodato & Rice 2004; Stamatellos & Whitworth 



20081 ; iForgan et al.ll2009h . the stellar motion is included automati- 



cally by treating the star as a central sink particle that is smoothed 
differently f rom the rest of the S PH particles. However, these simu- 
lations, as in lRice et alj ( l2003ar) , can exhibit a large initial accretion 
rate onto the central object. How the transfer of angular momentum 
from the disk material to the central object is handled could be quite 
important to the star's motion. The only global 3D grid-based hy- 
drodynamics scheme other than the lUHG code so far used for pub- 
lished simulations of GIs in protoplanetary disks around solar-type 
stars is the Eulerian spherical-grid code of Alan Boss. Althou gh he 
does not explicitly integrate the stellar equation of motion, IBossI 
( l2000 h allows the central protostar to wobble in response to the 
growth of nonaxisymmetry in the disk by repositioning the star to 
keep the system COM at the grid center. There is no guarantee with 
this scheme that the star's displacement is a true response to New- 
tonian reaction forces applied by the disk. The star's motion could 
simply represent an accumulation of numerical error in the loca- 
tion of the COM. In principle, then, this treatment is no better than 
what was done in the lUHG simulations. In order to explore the ef- 
fect of the star/disk interactio n, we have now im plemented the in- 
direct potential method (e.g.. iNelson et alJuOOOl) . This effectively 
puts our simulations into the accelerated reference frame of the star 
and therefore properly accounts for stellar motion while keeping 
the star at the center of our computational grid. 

In this article, we present the results of a simulation for a non- 
fragmenting disk using the indirect potential to treat stellar motion 
and compare it to results of an identical simulation with an arti- 
ficially fixed central star from Mejia et al. (2005). In SJS] we first 
present an overall qualitative comparison, followed by a detailed 
analysis of the stellar motion and of differences in the disk behav- 
iors. We compare our results to those obtained via other numerical 
methods and briefly discuss possible observational consequences 
in i|4] Finally, ^presents our main conclusions. 



2 METHODOLOGY 

2.1 The lUHG hydrodynamics code 

The two disk simulations in this paper are evolved using the 
standard lUHG 3D hydrodynamics code dPickett et alj 1200 3l ; 



iMeiia et alj|2005l) . This code solves the equations of hydrodynam- 
ics in conservative form on an evenly spaced Eulerian cylindrical 
grid (zo, (j), z) to second order in space and time. For these calcula- 
tions, the equation of state is that of an ideal gas with a ratio of spe- 
cific heats 7 = 5/3. The rotation axis is in the 2-direction, reflec- 
tion symmetry is assumed about the equatorial plane, and the cell 
widths in the vj and z-directions are the same. Shocks are medi- 
ated through the inclusion of von Neumann- Winkler artific ial bulk 
viscosity terms in the momentum and energy equations iPicketll 
Il995l) . and free outflow boundary conditions are assumed along the 
top, outer, and inner (central hole) edges of the grid. The (■co,(j),z)- 
directions will hereafter be referred to as the cylindrically radial, 
azimuthal, and vertical directions, respectively. When needed, the 
spherical radius from grid center is represented by r. Cooling is 
introduced into the energy equation as a local volumetric cooling 
rate A = e/icooi, where e is the internal energy density and icooi 
is a constant time scale on which the disk is assumed to be losing 
thermal energy due to radiation. 

The potential whose gradient yields the gravitational source 
terms in the momentum equation of the disk is made up of several 
components, namely, the disk component, the stellar component, 
and the indirect component. So 



<E>tot = $disk + "l-star + $1,, 



(1) 



$disk is calculated by first determining the disk boundary potential 
by multipole expansion of spherical harmonics up to / = \m\ — 10. 
Once the boundary conditions have been obtained, the density data 
is Fourier transformed in the (^-direction. Each of the Fourier com- 
ponents yields a 2D boundary value problem which can be cast 
into a block-tridiagonal matrix and solved using cyclic reduction 
(Tohline 1980). The solution in Fourier space is then transformed 
back into real space. "I>star is — GA/,/r, and <l>ind is described 
in the next section. For the comparison calculation in lMeiiaet al.l 
( 120050 ■ "IJind is zero. Although the star remains at the grid center in 
both simulations that we discuss, we refer to the lMeiia et al.l J2005h 
simulation without the indirect potential as the "fixed star" case and 
the simulation with the indirect potential as the "indirect case". 

In our earlier constant icooi simulations with 256 or 512 ra- 
dial grid cells that have a fixed central stellar potential, the disk's 
COM moved at most only a few radial cells away from the grid cen- 
ter over many disk rotations due to numerical errors. Nevertheless, 
even if fixing the star at the center of the inertial frame does not 
lead to pathological behavior, it also does not capture the full inter- 
action between the star and disk, including possible feedback be- 
tween stellar motion and growth of spiral modes (e.g.. I Adams et al.l 
il989 ; Shuetal. 1990) . 



2.2 Indirect potential 

To model the stellar motion explicitly, we use the indirect potential 
method (e.g., lNelson et alj200tt) . The grid is now considered to be 
the accelerated reference frame of the star, with the practical benefit 
that the star remains fixed at the grid center. The acceleration of 
the star-centered reference frame due to gravitational forces can be 
included into the fluid equation of motion as an a dditional grad ient 
of a potential, termed the indirect potential. As in lNelson et all the 
indirect potential at some point in the disk is given by 



•l-i^ 



Gl ^^vr', 



(2) 
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Figure 1. Midplane densities in logarithmic scale of the indirect (top) and fixed (bottom) simulations. The axes have units of AU and the time in the upper 
right is given in ORPs. The small white x in the indirect panels is the position of the system center of mass (disk plus star). The left panels show the saturation 
of the inidal growth in the burst phase. The right panels are representative of the asymptotic behavior. 



where the integration over the primed coordinates extends over the 
whole computational grid. In practice, we compute 



and then 






Gr- I ^^r' 



(3) 



(4) 



for each cell. Of course, this turns out to have much less computa- 
tional burden than computing the entire integral in equation lO for 
each cell. 



2.3 Initial model 



For th is comparison, we use the same initial disk as iMeiia et al.l 
1 120051) , except that we rescale the physical parameters so that 
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Matar = 1 M© and Mdisk = 0.14 Mq. The lMeii'a et al J simulation 
results are also appropriately rescaled in all comparisons presented 
here. The disk extends from 2.3 to 40 AU with a surface density 
E oc tj7~^''^. The computational gird has a central hole with a ra- 
dius of 1.6 AU in which the star resides, and the 40 AU initial disk 
outer radius is at radial zone 242. The new run with the indirect 
potential is started from the same axisymmetric disk after applying 
the same 0.01% amplitude random cell-to-cell density perturba- 
tion used by^ejia et al.. As in Meii a et al.. the constant icooi is set 
equal to 2 Outer Rotation Periods or ORPs. An OR? is defined as 
the initial rotation period, about 180 yr, at radial zone 200 (~ 33 
AU). The cylindrical grid in {■uj,4>,z) is initially 256x128x32. An 
additional 256 radial zones are added at about 5 ORPs to extend 
the grid to about kj = 85 AU once GIs cause the disk to expand ra- 
dially off the initial computational grid. The initial Q-distribution 
has a minimum of about 1 .5 at tj7 ~ 34 AU, and so the disk is only 
marginally stable at time i = 0. The indirect potential simulation 
is run for a total of about 19.6 ORPS or 3,500 yr. 

In the iMeiia et alj ( 120050 simulation, the fixed initial central 
"star" is a radially extended oblate mass distribution inside the cen- 
tral disk hole. When fitting the Mejia star by a point mass for the 
indirect case, an error of 0.1% is made in the star's mass. An ad- 
ditional error in the central force due to the nonsphericity of the 
Mejia star is about 0. 1 % at 11 AU and falls off quickly outside that 
radius. These discrepancies are two orders of magnitude less than 
the tens of percent differences reported for the GI behavior outside 
10 AU and should be unimportant. 



3 RESULTS 

3.1 Overall evolution 

3.1.1 Evolutionary phases 

Qualitatively, the overall outcomes are fairly similar. The disk sim- 
ulated with the indirect potential goes through the same four phases 
as the disk in the fixed star case iMeiia et al. 2005) . namely axisym- 
metric cooling, the onset of GIs in a burst, and a transition to a final 
quasi-steady asymptotic phase. The dividing times between these 
phases are also roughly the same. The onset of the burst phase oc- 
curs around 4.5 ORPs with the initial burst being predominately in 
discrete four to six-armed modes. The transition phase begins near 
7 ORPs and continues until the start of the quasi-steady asymptotic 
phase around 11 to 12 ORPs, after which heating by GI activity 
balances cooling on average throughout the Gl-active region and 
instability of a large number of interactin g spiral waves is sustained 
in an approximate steady state (see also .Gammiell200 ih . 



3.1.2 Final state 

The right-hand panels in Figure [T] offer a comparison of the final 
midplane densities for the disks. The structures in the indirect po- 
tential simulation appear to be a bit less sharply defined, probably 
owing to somewhat weaker GI amplitudes, as discussed in i]3.1.4| 
One can also see that the system COM, indicated by the white x 
in Figure [T] has remained within the central hole of the grid. Al- 
though this motion is not large from the disk's perspective, it may 
be of observational interest (see ^A.2\ . Additionally, as can be seen 
in Figure |2] the final azimuthally-averaged surface density profiles 
show very little difference, except that the indirect potential disk is 
slightly less extended radially. Assuming the surface density profile 
is a power law E oc r"'', a least squares fit for q between 15 and 



40 AU gives 2.33 and 2.31 for the indirect potential and fixed star 
cases, respectively. Also, the Toomre Q-values are similar in the 
asymptotic phase and average to about 1.2 to 1.3 over the 15 to 40 
AU region. 



3.1.3 Differences in the burst phase 

Despite the overall similarities, the GI structure for the indirect case 
exhibits some interesting differences. These can best be illustrated 
quantitatively by examining the overall amplitudes Am of sin(7n</f>) 
and cos{m(l>) terms in a Fourier decomposition of the density in 
the azimuthal directio n. We compute the Fourier components as in 
llmamura et alj ( 120001) such that 



^m — 



TT J pordrdz 



where 



and 



a-m = / pcos{m(f>)rdrdzd(f) 



hm = 



/ ps\n{m(t>) 



rdrdzdtj). 



(5) 



(6) 



(7) 



Here po is the axisymmetric component of the density, and the inte- 
grals extend over the computational grid. We typically sample the 
value of the components 100 times per ORP and, for some pur- 
poses, average them over a time interval, which is typically several 
ORPs in the asymptotic phase. 

As shown in Figure[T] there is a qualitative difference in dom- 
inant mode during the burst. For both cases, analyses of Am{t) 
plots reveal exponential growth in the linear regime at similar rates 
for m — 3 to 6 spiral waves centered near 30 AU, the region 
where Q is initially minimum. On the other hand, the m — 4 
and 5 waves play very different roles in the two cases. A4 grows 
first and is dominant for the fixed case; growth in A5 is substan- 
tially delayed. On the other hand, A5 and A^ dominate growth for 
the indirect case, with m = 4 remaining unimportant until ampli- 
tudes become nonlinear. All spiral waves grow from random un- 
correlated noise imposed initially in both cases that displaces the 
COM from the grid center. In the fixed case, where deviations of 
the star from the COM are not handled physically, an ordering in 
the small (~ 10~*) amplitudes of m = 4, 5, and 6 develops within 
a small fraction of an ORP such that m — 4 has a head start in 
its linear growth phase and reaches nonlinear amplitude about a 
third of an ORP before the other modes. On the other hand, in the 
indirect case, where the COM is treated properly, the early val- 
ues of A2 to Aq are more nearly the same, as they should be for 
white noise. Then m = 5, followed closely by m, = 6, dominates 
the linear growth. During the linear growth phase for both cases, 
Al ~ ^43^4 + AiA^ + ^5^46, as expected if m = 1 results from 
the nonlinear interaction of the rn = 3 to 6 modes and is not itself 
an independent unstable mode (e.g., Laughlin & Korchagin 199q)- 
As discussed in ^^^ the motion of the star with respect to the COM 
during the burst of the indirect case appears to be a response to this 
nonlinear interaction and should be sensitive in detail to the modes 
involved. Because our initial conditions are arbitrary and the subse- 
quent evolution is not substantively different, the qualitative differ- 
ence in dominant mode during linear growth does not undermine 
the general conclusions about disk evolution drawn from simula- 
tions with fixed stars in earlier lUHG papers. On the other hand, the 
significant displacement of the star during the burst and subsequent 
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phases is a potentially interesting and observable consequence (see 

E3. 



Fixed (Block) Indirect (Red) 



T = 0.15 
T = 19.57 
T = 19.57 



3.1.4 Differences in the asymptotic phase 

As shown in Figure[3] all the Am values are lower (about 0.05 to 0.2 
dex in Am, or on the order of tens of percent) for the indirect case 
during the asymptotic phase. The bigger differences appear to be 
in the higher-order Fourier components. Interestingly, although the 
indirect case allows for the possibility of SLING amplification, Ai 
is in fact measurably lower in the indirect case. The "error bars" 
on each of the points indicate the rms fluctuations over the sam- 
pling interval. It should be noted that the number of epochs used 
to compute the average Am and their fluctuations is different for 
each of the runs, namely 801 and 504 for the indirect and fixed star 
cases, respectively. This is due to sampling data for a fixed number 
of time steps, even though the time step length actually differed be- 
tween the calcualtions. Because the number of epochs sampled is 
large in both cases, we do not think this causes significant differ- 
ences in the average Am or rms fluctuations. 

As a measure of the total amplitude of non-axisymmetric den- 
sity structure, we can sum Am over m from 1 to 64. For the time 
period between 4.5 and 1 1 ORPs, representing the burst and transi- 
tion phases, the peak summed amplitude is 5.8 in the indirect case 
as opposed to 7.6 in the fixed star case. The average of the summed 
amplitudes over this time interval is also lower, at 2.1 for the in- 
direct case compared to 2.7 in the fixed star case. This trend con- 
tinues throughout the asymptotic phase with the mean value of the 
summed amplitudes being 2.2 in the indirect case compared to 2.7 
in the fixed case averaged over the interval from 12 to 19.5 ORPs. 
Figure [3] shows the full range of m values averaged over this time. 



3.1.5 Mass transport 

Figure |5] compares the disk mass transport rates during both the 
burst phase and the asymptotic phase, as determined by changes in 
the distribution of disk mass with respect to cylindrical radius vj 
measured from the star. Overall, we see that the addition of stel- 
lar motion does not affect the general outcome of GI activity for 
mass transport, even though it does produce a measurable weak- 
ening of the GI amplitudes. In both cases, typical radial inflow 
rates in the disks during the burst and transition phase are ~ 3 to 
7 X 10~^MQyr~^ over 15 to 35 AU. The inflow is harder to char- 
acterize for the asymptotic phase but crudely averages to something 
like '^ 7 X 10~'^M (7)yr~^ over 10 to 35 or 40 AU in both cases. 
iMeiia et al.l ( l2005h have already demonstrated how much fluctua- 
tion there is in mass inflow due to GIs (see the light grey curves in 
their Figure 3). 



3.2 Stellar motion 

With the indirect potential, we can infer the physical motion of the 
star by determining the position of the system COM on our grid 
as a function of time. Figure |6] shows the COM's equatorial plane 
trajectory in the star's reference frame and the distance ■oucomit) 
between the star and COM as a function of time for the duration of 
the simulation. In Figure|6] the black depicts the COM motion from 
to 12 ORPS, i.e., through the burst and transition phases, while 
the red shows the asymptotic phase for t >12 ORPs. The star does 
not become significantly displaced from the system COM until the 
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Figure 2. Surface densities in logarithmic scale of the indirect (red) and 
fixed (black) simulations. 



burst commences at 4.5 AU and achieves its largest excursion from 
the origin, ~ 0.25 AU, near the end of the burst at about 7 ORPs. 

As discussed in Section |3.1.3l there does not appear to be sus- 
tained growth of an independent m = 1 mode through a feedback 
loop. The initial displacement of the star during the burst appears 
to be due to the nonlinear interaction of odd and even spiral waves 
with 771 = 3 to 6. As can be seen in Figure|4]there is a clear corre- 
lation between the exponential growth of the Am and the displace- 
ment of the COM. The COM begins with a very small displace- 
ment due to the random perturbation in the initial disk. The COM 
remains at this small radius until about a half ORP after the onset 
of exponential growth in the Am&. The time lag is probably due to 
the need for a wave to propagate to the ed ge of the disk be fore the 
COM displaces (see step 3 in Section Ic of lShu et alj ( 119901) ). Once 
the disk enters the asymptotic phase, the COM excursion decreases 
in amplitude but remains non-negligible. From 12 to 19.5 ORPs the 
mean distance from the star is 0.09 AU, while the average distance 
of the COM from the star is about 0. 1 1 AU over the entire time 
period between 4 and 19.5 ORPs. 

To determine the periodicity of the COM motion, and hence, 
by reflex, of the star, we construct Lomb-Scargle periodograms 
JLombll 197^ : ISca rgle' 1982) for the equatorial plane COM motion 
in u7com and in Cartesian x and y coordinates. Figure |7] shows the 
power spectrum for x measured between 12 and 19.5 ORPs; the 
y periodogram (not shown) is similar. If we assume that the most 
powerful signal is due to the star's primary orbit motion, we get 
a period of 0.80 ORP (about 140 yr) for both x and y. The peri- 
odogram for tUcom gives consistent results but suggests a slightly 
shorter period of 0.72 ORP (about 130 yr). One expects the x 
and y periods to agree, while the zucom period should be shifted 
somewhat due to apsidal precession, as observed. Several statisti- 
cal tests were done to verify the significance of the signals, includ- 
ing changes in the sampling rate of the COM position (see Michael 
2010, Ph.D. dissertation). The signal is robust, and one also ob- 
serves a similar period just by counting peaks in the zucom{t) plot 
over the same time interval. As in Rice et al. (2003a), there are sig- 
nificant higher-order motions superposed on the basic COM wob- 
ble. 

During the burst phase, the period of the COM motion ap- 
pears to be somewhat longer, closer to 1 ORP, in rough agreement 
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Mass Flow Comparison, 4.39--10.96 ORP 
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Figure 3. Comparison of indirect (red) and fixed (black) Am values for 
eacli m averaged from 12 to 19.5 ORPs. The error bars represent the RMS 
of the fluctuations over this time period. 



Mass Flow Comparison, 12.09--19.58 ORP 





R (AU) 



Figure 5. Mass transport rates in solar masses per year plotted during the 
burst (top) and asymtotic (bottom) phases. The horizontal axis is given in 
AU and the red curve is the indirect simulation, the black curve is the fixed 
simulation and the dashed line represents no mass flow. 



Figure 4. Am (t) and t<7com (i) for the initial phase of the indirect simula- 
tion. The vertical Am scale is on the left; the COM scale is on the right. 



with the pattern period of the predominant five-armed mode (see 
the top left panel of Figure [Hi. Although the burst is dominated by 
the m = 5 pattern, both odd and even modes with similar pattern 
periods grow together. The resulting m = 1 distortion of the disk 
due to the nonlinear interaction of these modes appears to kick the 
star into orbit around the system COM with a period similar to the 
pattern period of the growing spiral waves. This is consistent with 
the following simple estimate. The star's largest excursion from the 
COM is 0.25 AU, while corotation for the bursting spiral waves is 
at about 30 AU. So the star's motion implies a net mass asymme- 
try in the disk of about 0.25Afstar/30 ^ O.OO8A/0 (about eight 
Jupiter masses) at about 30 AU on the other side of the COM. The 
mass in the disk between the Lindblad resonances of the dominant 
TO = 5 mode during the burst is approximately O.OSAf©. So just 
a fraction of the mass (a bit more than one spiral arm's worth) in- 
volved in the interacting odd (m — 3 and 5) and even (to — 4 and 
6) waves is sufficient to displace the star by the amount observed. 

The significant motions of the star in the inertial frame show 
that an appreciable amount of angular momentum is transfered be- 



tween the disk and the star. Figure[8]plots the star's angular momen- 
tum produced by the star/disk interactions, as calculated by taking 
a five point numerical time derivative of the system COM position 
vector to get a velocity vector. The COM position is determined 
by a numerical integral over the grid, and the velocities are then 
determined by numerical derivative, so numerical errors are ampli- 
fied and the result is fairly noisy. We have tried to estimate grav- 
itational torques by doing a numerical derivative of the curve in 
Figure[8l but the results have even more noise, and we are not sure 
how much of the fluctuation is numerical or due to the gravitoturbu- 
lence in the disk. The average torque during the asymptotic phase 
is a few x 10"'** ergs, with a suggestion that instantaneous torques 
can be much larger or smaller. 



3.3 Analysis of the disk structure 

In Figure [3] we compared the GI amplitudes in various Fourier 
components for the two simulations during the asymptotic phase. 
Another important characteristic of m-ar med spiral waves i s the 
coherence of their patterns. As explained in lMeiia et al J ( 120050 , we 
can search for patterns that are coherent over a range of radii by cre- 
ating power spectra at all radii for cos(</)m), where (pm is the phase 
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Figure 6. System COM motion. The top plot shows the COM trajectory 
in the equatorial plane, and both axes are labelled in AU. Black represents 
motion prior to 12 ORPs, and red after The bottom plot shows the system 
COM radial excursion for the duration of the simulation, radial excursion 
labelled in AU and time in ORPs. 
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Figure 8. Plot of the stellar angular momentum as a function of time as 
determined from the motion of the system COM. 



As pointed out in lMeiia et al.l d2005h . coherent patterns tend to be 
confined between their Lindblad resonances. 



of the m-th Fourier component of the azimuthal density structure. 
We restrict the analysis to the equatorial plane. Figure [9] compares 
periodograms for m, = 1, 2, and 3 for the two simulations during 
the asymptotic phase. Within each diagram, the relative power, as 
denoted by the color scale, is a measure of the coherence of the pat- 
terns detected at each radius. The apparent difference in color scale 
between plots for the two simulations is an artifact of somewhat 
different sampling rates with time and has no significance. What 
matters are the strong vertical strips, which indicate periods that 
are present over a range of radii, i.e., spatially coherent wave pat- 
terns. Curves representing the approximate locations of the outer 
Lindblad resonance (OLR), corotation (CR), and the inner Lind- 
blad resonance (ILR) for a given pattern period Ppat are shown. 
The curves are a bit ragged due to the use of azimuthal averages 
and/or numerical derivatives to compute these radii numerically. 



3.3.1 One-armed structure 

In Figure 121 m = 1 does not exhibit any stronger radially coher- 
ent patterns in the indirect simulation. In fact, except perhaps for 
some slow-moving distortion of the outer disk between 50 and 60 
AU, where there is little mass, neither simulation shows any vertical 
stripes of power in m = L There are no coherent global one-armed 
spiral pattens in the dis ks. Because SLING is most likely to am- 
plify one-armed modes JAdams et al]|l989l: IShu et al.lll990l) , this 
suggests that sustained SLING amplification is not present during 
the asym ptotic phase of eith er si mulation. Using the methods de- 
scribed in lBolev et aT] ( l2006h and lCai et aljjiooi) . FigurefTOlcom- 
pares the gravitational torque of the outer disk on the inner disk for 
several m-values at all -hj during the asymptotic phase. For both 
simulations, these disk internal gravitational stresses are measured 
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Table 1. Pattern Speeds of Prominent Periodogram Features 



Signal Measured 


Pattern Speeds (1/ORP) 


COM z-motion 


0.55 


0.80 


1.05 1.50 


COM j;-motion 


0.51 


0.80 


1.02 1.50 


m = 1 


0.65 


0.85 


— 1.50 


m = 2 


0.50 


0.80 


— 1.50 


m = 3 


0.59 


0.80 


1.05 1.50 



using cylindrical radii centered on the star. Figure [TO] shows negli- 
gible contribution to the torque from the m, = 1 component in both 
simulations. 
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3.3.2 Disk/star interaction 

Although the motion of the star does not induce sustained growth 
of coherent m = 1 structure, the star clearly has a strong interac- 
tion with the disk. Table [T] compares the strongest periods in the 
COM X and y periodograms and the strongest stripes of power in 
the disk periodograms. The periodicity in the x and y COM mo- 
tion corresponds to maxima in m = 1,2, and 3 periodograms 
at pattern speeds near 0.5, 0.8, and 1.50/ORP. The most striking 
features in the m — 2 and 3 periodograms are those at pattern 
speeds of 1.50/ORP and 1.05/ORP, respectively. Figure [TO] shows 
the effect of this star/disk interaction on the disk torque. Clearly, 
the disk torque is decreased in the indirect simulation when com- 
pared to the fixed simulation, especially near 24 AU, which cor- 
responds to the corotation radius of a pattern with a pattern speed 
of about 1.50/ORP. Although the interpretation is unclear, we also 
note that the torque maximum at 19 AU roughly corresponds to the 
inner Linblad resonance of a coherent m — 3 pattern with a pattern 
speed of near 1.5/ORP. Moreover, the decrease in disk torque cen- 
tered near 24 AU is comparable to the torque on the star. The time 
averaged torque from 12 to 19.5 ORPs on the star is about 5 x 10^^ 
ergs, while the difference in disk torques averaged over the same 
time interval and averaged from 10 to 60 AU is about 9 x 10^* 
ergs. This is essentially agreement within the accuracy to which 
these averaged torques can be measured. Apparently, this amount 
of torque, which is confined to the disk in the fixed star case, is 
effectively transferred to the star when the star is allowed to move. 
Even though these differences in the disk structure introduced 
by the star's motion are measurable in the torque profiles and pe- 
riodograms, the overall mass transport in the disk is largely unaf- 
fected as can be seen in the mass transport rates of Figure [5] An- 
other metric of mass transport is the effective gravitational a, which 
is plotted for both the fixed and in direct simulatio ns in Figure [TT] 
FigurefTTjalso shows the predicted .Gammig ( 1200 ih values for fully 
self-gravitating and non-self-gravitating disks for tcooi = 2 ORP 
when the disks are in a local balance between radiative cooling and 
heating by GI activity. Although the indirect potential a is some- 
what smaller than the fixed-star q, it follows the overall trend of 
the fixed-star a. The decrease in a is consistent with the decrease 
in torque, b ecasue these effect ive a-values are computed from the 
torques (see Bolev et alJl200q) . Notice that the inclusion of stellar 
motion brings the measured a so mewhat closer to values predicted 
by the local balance assumption JGammiell200lh . which are also 
shown on the plot. The agreement between the a-values from in- 
direct simulation and the local treatment are closest in the region 
around the corotation of the strongest m — 2 and 3 patterns in Fig- 



Figure 10. Torque profiles averaged from 12 to 19.5 ORPs for the indirect 
(red) and fixed (black) simulations. 




Figure 11. Effective Shakura-Sunyaev Q-values computed for the fixed star 
(black curve) and indirect (red curve) simulations averaged over the asymp- 
totic phase from 12.5 to 19 ORPS. Show n for comparison are curves pre- 
dicted by equation (13) in lGammid 1200 ih with tcool = 2 ORP. The upper 
curve assumes a vetical structure dominated by disk self-gravity; the lower 
curve assumed the vertical gravity is due to the star alone. For Q fa 1.3, we 
expect the disk to be significatly self-gravitating. 



ure[9] It was argued bv lBalbus & Papaloizoul <1999f) that agreement 
between global and local energetic behavior of GIs was most likely 
to occur near corotation. 



4 DISCUSSION 

4.1 The effect of star/disk interactions 

We have found that the inclusion of stellar motion in our 3D ra- 
diative hydrodynamic simulations has a measurable, but not drastic 
effect on the evolution of a gravitationally unstable protoplanetary 
disk. Gross features such as surface density profile. Am values, and 
mass transport rates remain unchanged within tens of percent. On 
the other hand, the motion of the central star is not negligible, and 
there is a relationship between changes in the disk self-torque and 
the star/disk interaction. 
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Figure 9. The m = 1,2, and 3 phase periodograms measured using data from 12 to 19.5 ORPs. The color scale represents relative power with purple being 
the smallest and red/white being the largest. In this plot, relative power is the important measure; the absolute numbers are in arbitrary units. The fixed star 
case is on the left, and the indirect case on the light. 
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During the burst phase, the freedom of the star to move qual- 
itatively changes the spiral modes that dominate at the initiation of 
GI activity, and the stellar motion reflects the orbit period at the 
CR of the growing modes. During the asymptotic phase, there is 
also a strong correlation between the quasi-periodic motions of the 
star and periodicities evident in the m = 1 — 3 structure, as seen by 
comparing Figures|7]and|9]with Table[T] Despite the complex inter- 
action between the disk and star that gives rise to this correlation, 
measures of Gl-activity in the disk during the asymptotic phase 
only tend to be reduced by at most about tens of percent. Overall, 
it appears that the net effect of the star/disk interaction is to intro- 
duce a new degree of freedom to the system which allows some of 
the torque that previously went into disk self-torque to be applied 
instead to the star/disk interaction. This, in turn, slightly lowers the 
Am values, disk self-torques, and resultant a. The star/disk torque 
interaction transfers a significant amount of angular momentum 
into orbital motion of the star. Further study is required to deter- 
mine the sensitivity of this interaction to various simulation condi- 
tions, such as cooling time (and thereby GI strength), initial surface 
density profile, and disk-to-star mass ratio. 

Figure [To] which divides the contribution to the disk torque 
into its various Fourier components during the asymptotic phase, 
shows negligible contribution from m = 1 or one-armed distur- 
bances, which would be expected to dominate if SLING were ef- 
fective in our disk. Analysis of the Am{t) plots during the burst 
(not illustrated here) also shows that m = 1 only grows as a result 
of the nonlinear interaction of spiral modes with higher m as the GI 
activity begins. There are two possible reasons for failure to detect 
evidence for a one-armed instability: 1) Although a dynamic pro- 
cess, the growth times for one-armed instabilities, whether driven 
by SLING or not, tend to be several to many outer rotation periods 
( lAdams et al.ll 19891 : IShu et alJI 19901) , because one-armed instabili- 
ties involve a reflection from the outer boundary of the disk. Growth 
times become especially long for disk-to-star mass rati os of order 
0.1, a s illustrated by the low growth rates in Figure 4a of lNoh et al.l 
( 11992 ). Our simulations probably do not last long enough for one- 
armed instabilities to grow. 2) SLING amplification in particular 
prob ably does not occur because our disk-to-total-mass ratio is too 
low JShu et alJll990l : lNoh et alj|l992h . 



4.2 Observable consequences of stellar motion 

The stellar displacement from the system COM due to star/disk in- 
teractions is as large as 0.25 AU during the burst and averages 0.09 
AU during the asymptotic phase in our simulation with the indirect 
potential. S uperficially , this se ems strikingly at odds with results 
reported by iRice et alj ( l2003ah . who found a maximum radial ex- 
cursion by the central star of ~ 5 x 10"'' AU . 

While it is true that the numerical methods ar e fairly differ- 
ent, in that we use a grid-based code and lRice et alj use SPH, most 
of the discrepancy can probab ly be attrib uted to several other key 
differences: 1.) Cooling time. i Rice et al.l use a local cooling time 
tied to the local orbital speed. Specifically, they present simula- 
tions for fcooi = 5f2~^ and tcooi = 10S1~^ or, equivale ntly. 
tcooi/Prot = about 0.8 and 1.6, respectively. The [RiceetalJ pre- 
scription has the effect of shortening the cooling time with decreas- 
ing radius, whereas the cooling time in our simulation is constant 
across the disk and is chose n to be twi ce the Prot at 33 AU. 2.) 
Disk Structure. Although the lRice et alj disk has a similar disk- to- 
star mass ratio, 0. 1 compared with our 0. 14, their inner disk has a 
much higher fraction of the mass, because their disk surface den- 
sity distribution is much steeper, roughly E ~ zu~'''^ compared 



with our E ~ c<7~^". 3.) Stellar Accretion. While very little mass 
moves across the cylindrical inner boundary of our grid, the SPH 
calculation allows its particles to approach the star, and they are 
incorporated into it once they enter its smoothing radius. 

Because of the short cooling times in the innermost parts of 
the lRice et akl disk simulations, the inner dis ks are as violently GI 
unstable as the ou ter disks. As discussed bv lMeiia et alj 1 120050 in 
comparisons with iLodato & Ricd j2004h . major consequences of 
using a constant tcooi/Piot prescription rather than a constant tcooi 
prescription in a moderate mass disk are that there is no burst, be- 
cause the whole disk goes unstable at once, and that the GIs are not 
dominated by a few individual global modes. As our analysis has 
shown, it is the interaction of the star with a few discrete modes 
in the outer disk in both the burst and the asympt otic phases that 
causes the large displacements. Also, in lRice et all the disks expe- 
rience a large amount of mass accretion. In fact, the inner 2 AU of 
their disks is rapidly accreted onto the central star, and they assume 
that accreted angular momentum goes into rotation, not orbital mo- 
tion. In our simulation, the inner disk remains stable, and there is 
essentially no accretion. The outcome of accretion onto the star, if 
it occurs in a simulation, could vary significantly depending on as- 
sumptions about how to treat the linear and angular momentum of 
any material added to the star in an asymmetrical ma nner. 

As with many other aspects of GI behavior (e.g., Pickett et al.l 
1 19981 I2OOOI l2003l: iMeiia e"tai]|2005l : Isolev et alj|2006l) . displace- 
ment of the star is evidently quite sensitive to assumptions about 
disk structure and, especially, disk thermodynamics. What our in- 
direct simulation shows is that, under conditions where GIs initiate 
in a burst and tend to be dominated by a few discrete modes in the 
outer disk, the stellar excursions about the COM can be quite large. 
In our simulation, the periods of the induced stellar orbital mo- 
tion correspond to the CRs of the dominant modes, or about 130 
and 170 years during the asymptotic and burst phases, respectively. 
The amplitudes of the displacements correspond to disk asymme- 
tries with masses on the order of several to many Jupiter masses. In 
the asymptotic phase of our simulation, the astrometric signature 
would be on the order of 16 /las/yr if the system were at a dis- 
tance of 100 pc, and the maximum radial velocity would be about 
30 m/sec. Stellar motion of this type could be difficult to distin- 
guish from that caused by one or more super-Jupiter protoplanets 
orbiting at lO's AU. 

Our simulation with the indirect potential demonstrates that 
discrete global GI spiral modes can, in principle, drive large stellar 
motion. However, whether this effect will be manifest in real disks 
may depend in part on whether accretion through the inner disk 
onto the star happens in an asymmetric fashion that counteracts 
the GI torques or whether stellar orbital angular momentum can be 
effectively removed by a jet or wind. Orbital motion of the star and 
of the innermost disk due to GIs might cause a precessing or helical 
jet/wind structure. These would be interesting questions for future 
research. 



5 CONCLUSIONS 

The two major goals of this paper were: 1) to determine how much 
a proper treatment of stellar motion affects qualitative and quan- 
titative conclusions made in previous studies by our group and 2) 
to decide whether stellar motion due to star/disk interaction could 
lead to observable consequences. We find that, although the simu- 
lations do show a weakening of disk self-torques and nonaxisym- 
metric structure by at most about tens of percent, simulations with 
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and without proper treatment stellar motion are overall very similar. 
They exhibit the same evolutionary phases, and they both settle into 
similar quasi-steady asymptotic states where cooling is balanced by 
GI heating. The main qualitative difference is a shift of dominance 
from an even mode (four-armed) to an odd mode (five-armed) dur- 
ing the burst phase when the star's motion is properly treated. This 
would only matter in a case where the initial conditions at the time 
of onset for GIs were well understood for a particula r application. 
We conclude that lUHG s tudies using a fixed star (e.g. JPickett et al.l 
20031: iMeiia et al.ll2005l; ICai et al.ll2006l : iBolev et aljr2006l l2007l : 



Bolev & Durisenll2008l : ICai et alj |2008') have produced valid and 



reliable results regarding the overall behavior and evolution of GI- 
active disks. 

When we account for stellar motion through an indirect po- 
tential, the motion of the star appears to be due to the nonlinear in- 
teraction of discrete odd and even global spiral modes, both during 
the burst and in the asymptotic phase. These interactions involve 
substantial disk mass at radii of about 30 AU (burst) and 24 AU 
(asymptotic phase) with the result that both astrometric and radial 
velocity disturbances are fairly large. For the asymptotic phase, we 
estimate about 16 /las/yr (at 100 pc) and 30 m/s, respectively. The 
time signatures can involve multiple periodicities, but it would be 
difficult to distinguish this from a stellar response to the presence of 
multiple massive planets or protoplanets. Our results are in contrast 
to the m uch more modest stellar motions reported by iRice et al.l 
( l2003an . under conditions where the cooling times decrease inward 
in the disk, there is no burst, GIs are not dominated by a few dis- 
crete global modes, and there is significant accretion from the inner 
disk. This means that stellar motion induced by GIs is sensitive 
to various aspects of disk structure and physics. Further study is 
clearly warranted, especially if outbursts of GIs are common during 
the early phases of disk forma tion when disks are large and mas- 
sive (e.g. J Laughlin & Bodenhe ime3l 1994IVorobvov & Basij|2006l : 
IStamatellos & Whitworthii2009l : lBolevll2009l) . 
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